function [b,a] = fir_win(N,Wn,varargin)
%FIR1   FIR filter design using the window method.
%   B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
%   and returns the filter coefficients in length N+1 vector B.
%   The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 
%   corresponding to half the sample rate.  The filter B is real and
%   has linear phase, i.e., even symmetric coefficients obeying B(k) =
%   B(N+2-k), k = 1,2,...,N+1.
%
%   If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
%   order N bandpass filter with passband  W1 < W < W2.
%   B = FIR1(N,Wn,'high') designs a highpass filter.
%   B = FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
%
%   If Wn is a multi-element vector, 
%          Wn = [W1 W2 W3 W4 W5 ... WN],
%   FIR1 returns an order N multiband filter with bands
%    0 < W < W1, W1 < W < W2, ..., WN < W < 1.
%   B = FIR1(N,Wn,'DC-1') makes the first band a passband.
%   B = FIR1(N,Wn,'DC-0') makes the first band a stopband.
%
%   For filters with a passband near Fs/2, e.g., highpass
%   and bandstop filters, N must be even.
%   
%   By default FIR1 uses a Hamming window.  Other available windows,
%   including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin
%   can be specified with an optional trailing argument.  For example,
%   B = FIR1(N,Wn,kaiser(N+1,4)) uses a Kaiser window with beta=4.
%   B = FIR1(N,Wn,'high',chebwin(N+1,R)) uses a Chebyshev window.
%
%   By default, the filter is scaled so the center of the first pass band 
%   has magnitude exactly one after windowing. Use a trailing 'noscale' 
%   argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'), 
%   B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale').
%
%   See also KAISERORD, FIRCLS1, FIR2, FIRLS, FIRCLS, CREMEZ,
%            REMEZ, FREQZ, FILTER.

%   FIR1 is an M-file implementation of program 5.2 in the IEEE
%   Programs for Digital Signal Processing tape. 

%   Author(s): L. Shure
%              L. Shure, 4-5-90, revised
%              T. Krauss, 3-5-96, revised
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.20 $  $Date: 1997/12/02 18:36:03 $

%   Reference(s):
%     [1] "Programs for Digital Signal Processing", IEEE Press
%         John Wiley & Sons, 1979, pg. 5.2-1.

nargchk(2,5,nargin);

% Up to 3 optional input arguments, always in this order:
%   1 - Filter type flag, can be 'high', 'stop', '', 'DC-0', or 'DC-1'
%   2 - Window vector
%   3 - 'noscale' flag

% default optional parameter values:
Ftype = '';
Wind = [];
Scale = 'scale';

switch length(varargin)
case 1
    if isstr(varargin{1})&(length(varargin{1})>0)
        s = upper(varargin{1});
        switch upper(s)
        case {'SCALE','NOSCALE'}
            Scale = s;
        otherwise
            Ftype = s;
        end
    else
        Wind = varargin{1};
    end
case 2
    if isstr(varargin{1})
        Ftype = varargin{1};
    else
        Wind = varargin{1};
    end
    if isstr(varargin{2})
        Scale = varargin{2};
    else
        Wind = varargin{2};
    end
case 3
    Ftype = varargin{1};
    Wind = varargin{2};
    Scale = varargin{3};
end

switch upper(Scale)
case 'NOSCALE'
    SCALING = 0;
case 'SCALE'
    SCALING = 1;
otherwise
    error('Scaling option must be ''noscale'' or ''scale''.')
end

if isempty(N) | ~isnumeric(N) | ~isreal(N) | N~=round(N) | N<=0
    error('N must be a real, positive integer.')
end

Ftype = upper(Ftype);
if ~strncmp(Ftype,'HIGH',1) & ~strncmp(Ftype,'STOP',1) & ...
   ~strncmp(Ftype,'DC-0',4) & ~strncmp(Ftype,'DC-1',4) & ...
   ~isempty(Ftype)
    error('Filter type must be ''high'',''stop'',''DC-0'', or ''DC-1''.')
end

nw = length(Wind);

nbands = length(Wn) + 1;
if (nbands > 2) & isempty(Ftype)
    Ftype = 'DC-0';  % make sure default 3 band filter is bandpass
end
First_Band = isempty(findstr('DC-0',Ftype)) & isempty(findstr('HIGH',Ftype));
mags = rem( First_Band + (0:nbands-1), 2);

L = N + 1;
odd = rem(L, 2);
if (mags(nbands) & ~odd)
      disp('For highpass and bandstop filters, order must be even.')
      disp('Order is being increased by 1.')
      N = N + 1;  L = L + 1;
      odd = 1;
end
if nw ~= 0 & nw ~= L
   error('The window length must be the same as the filter length.')
end
if nw == 0   % replace the following with the default window of your choice.
   Wind = hamming(L);
end
%
% to use Kaiser window, beta must be supplied
% att = 60; % dB of attenuation desired in sidelobe
% beta = 0.1102*(att-8.7);
% wind = kaiser(L,beta);

if  any( Wn<0 | Wn>1 )
   error('Frequencies must fall in range between 0 and 1.')
end
if  any(diff(Wn)<0)
   error('Frequencies must be increasing')
end

Wn = Wn(:)';
ff = [0,Wn(1:nbands-1); Wn(1:nbands-1),1];
mags = [mags(:)'; mags(:)'];
hh = firls(L-1,ff(:),mags(:));

b = hh.*Wind(:)'; 
a = 1;

if SCALING
    if First_Band
        b = b / sum(b);  % unity gain at DC
    else
        if ff(4)==1
            % unity gain at Fs/2
            f0 = 1;
        else
            % unity gain at center of first passband
            f0 = mean(ff(3:4));
        end
        b = b / abs( exp(-j*2*pi*(0:L-1)*(f0/2))*(b.') );
    end
end

function w = chebwin(n_est, r)
%CHEBWIN Chebyshev window.
%    W = CHEBWIN(N,R) returns the N-point Chebyshev window 
%        with R decibels of ripple.

%   Author: James Montanaro
%   Reference: E. Brigham, "The Fast Fourier Transform and its Applications" 

[n,w,trivalwin] = check_order(n_est);
if trivalwin, return, end;

w = chebwinx(n,r);


function w = hamming(n_est,sflag)
%HAMMING Hamming window.
%   W = HAMMING(N) returns the N-point symmetric Hamming window 
%       in a column vector. 
%   W = HAMMING(N,SFLAG) generates the N-point Hamming window 
%       using SFLAG window sampling. SFLAG may be either 'symmetric' 
%       or 'periodic'. By default, 'symmetric' window sampling is used. 


error(nargchk(1,2,nargin));

[n,w,trivalwin] = check_order(n_est);
if trivalwin, return, end;

% Set sflag to default if it's not already set:
if nargin == 1,
   sflag = 'symmetric';
end

switch lower(sflag),
case 'periodic'
   w = sym_hamming(n+1);
   w = w(1:end-1);
case 'symmetric'
   w = sym_hamming(n);
otherwise
	error('Sampling must be either ''symmetric'' or ''periodic''');
end

function w = sym_hamming(n)
w = .54 - .46*cos(2*pi*(0:n-1)'/(n-1));

function [n_out, w, trivalwin] = check_order(n_in)
%CHECK_ORDER Checks the order passed to the window functions.
% [N,W,TRIVALWIN] = CHECK_ORDER(N_ESTIMATE) will round N_ESTIMATE to the
% nearest integer if it is not alreay an integer. In special cases (N is [],
% 0, or 1), TRIVALWIN will be set to flag that W has been modified.

%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.5 $  $Date: 1997/11/26 20:14:55 $

w = [];
trivalwin = 0;

% Special case of negative orders:
if n_in < 0,
   error('Order cannot be less than zero.');
end

% Check if order is already an integer or empty
% If not, round to nearest integer.
if isempty(n_in) | n_in == floor(n_in),
   n_out = n_in;
else
   n_out = round(n_in);
   warning('Rounding order to nearest integer.');
end

% Special cases:
if isempty(n_out) | n_out == 0,
   w = zeros(0,1);               % Empty matrix: 0-by-1
   trivalwin = 1; 
elseif n_out == 1,
   w = 1;
   trivalwin = 1;   
end

function [h,a]=firls(N,F,M,W,ftype);
% FIRLS Linear-phase FIR filter design using least-squares error minimization.
%   B=FIRLS(N,F,A) returns a length N+1 linear phase (real, symmetric
%   coefficients) FIR filter which has the best approximation to the 
%   desired frequency response described by F and A in the least squares 
%   sense. F is a vector of frequency band edges in pairs, in ascending 
%   order between 0 and 1. 1 corresponds to the Nyquist frequency or half
%   the sampling frequency. A is a real vector the same size as F 
%   which specifies the desired amplitude of the frequency response of the
%   resultant filter B. The desired response is the line connecting the
%   points (F(k),A(k)) and (F(k+1),A(k+1)) for odd k; FIRLS treats the 
%   bands between F(k+1) and F(k+2) for odd k as "transition bands" or 
%   "don't care" regions. Thus the desired amplitude is piecewise linear
%   with transition bands.  The integrated squared error is minimized.
%
%   B=FIRLS(N,F,A,W) uses the weights in W to weight the error. W has one
%   entry per band (so it is half the length of F and A) which tells
%   FIRLS how much emphasis to put on minimizing the integral squared error
%   in each band relative to the other bands.
%   
%   B=FIRLS(N,F,A,'Hilbert') and B=FIRLS(N,F,A,W,'Hilbert') design filters 
%   that have odd symmetry, that is, B(k) = -B(N+2-k) for k = 1, ..., N+1. 
%   A special case is a Hilbert transformer which has an approx. amplitude
%   of 1 across the entire band, e.g. B=FIRLS(30,[.1 .9],[1 1],'Hilbert'). 
%
%   B=FIRLS(N,F,A,'differentiator') and B=FIRLS(N,F,A,W,'differentiator')
%   also design filters with odd symmetry, but with a special weighting
%   scheme for non-zero amplitude bands. The weight is assumed to be equal 
%   to the inverse of frequency, squared, times the weight W. Thus the 
%   filter has a much better fit at low frequency than at high frequency. 
%   This designs FIR differentiators.
%
%   See also REMEZ, FIR1, FIR2, FREQZ and FILTER.

%   Example of a length 29 low-pass filter:
%   	h=firls(30,[0 .1 .2 .5]*2,[1 1 0 0]);
%   Example of a low-pass differentiator: 
%   	h=firls(N,[0 10 20 40]/50,[0 .2 0 0],'differentiator');
%
%       Author(s): T. Krauss
%   History: 10-18-91, original version
%            3-30-93, updated
%            9-1-95, optimize adjacent band case
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.25 $  $Date: 1997/12/02 19:16:50 $

% check number of arguments, set up defaults.
nargchk(3,5,nargin);
if (nargin==3),
    W = ones(length(F)/2,1);
    ftype = '';
end
if (nargin==4),
    if isstr(W), 
        ftype = W; W = ones(length(F)/2,1);
    else
        ftype = '';
    end
end	
if (nargin==5),
    if isempty(W), 
        W = ones(length(F)/2,1);
    end
end
if isempty(ftype)
    ftype = 0;  differ = 0;
else
    ftype = lower(ftype);
    if strcmp(ftype,'h')|strcmp(ftype,'hilbert')
        ftype = 1;  differ = 0;
    elseif strcmp(ftype,'d')|strcmp(ftype,'differentiator')
        ftype = 1;  differ = 1;
    else
        error('Requires symmetry to be ''Hilbert'' or ''differentiator''.')
    end
end

N = N+1;                   % filter length
F=F(:)/2;  M=M(:);  W=sqrt(W(:));  % make these guys columns
dF = diff(F);
if (rem(length(F),2)~=0)
    error('F must have even length');
end;
if (length(F)~=length(M))
    error('F and A must be equal lengths');
end;
if (length(F)~=length(W)*2)
    error('There should be one weight per band.');
end;
if any(dF<0),
    error('Frequencies in F must be nondecreasing.')
end
if (max(F)>1)|(min(F)<0)
    error('Frequencies in F must be in range [0,1].')
end
if all(dF(2:2:length(dF)-1)==0)
    fullband = 1;
else
    fullband = 0;
end
if all((W-W(1))==0)
    constant_weights = 1;
else
    constant_weights = 0;
end

L=(N-1)/2;

Nodd = rem(N,2);

if (ftype == 0),  % Type I and Type II linear phase FIR 
                  % basis vectors are cos(2*pi*m*f) (see m below)
    if ~Nodd
        m=(0:L)+.5;   % type II
    else
        m=(0:L);      % type I
    end
    k=m';
    need_matrix = (~fullband) | (~constant_weights); 
    if need_matrix
        I1=k(:,ones(size(m)))+m(ones(size(k)),:);    % entries are m + k
        I2=k(:,ones(size(m)))-m(ones(size(k)),:);    % entries are m - k
        G=zeros(size(I1));
    end

    if Nodd
        k=k(2:length(k));
        b0=0;       %  first entry must be handled separately (where k(1)=0)
    end;
    b=zeros(size(k));
    for s=1:2:length(F),
        m=(M(s+1)-M(s))/(F(s+1)-F(s));    %  slope
        b1=M(s)-m*F(s);                   %  y-intercept
        if Nodd
            b0 = b0 + (b1*(F(s+1)-F(s)) + m/2*(F(s+1)*F(s+1)-F(s)*F(s)))...
                         * abs(W((s+1)/2)^2) ;
        end
        b = b+(m/(4*pi*pi)*(cos(2*pi*k*F(s+1))-cos(2*pi*k*F(s)))./(k.*k))...
                     * abs(W((s+1)/2)^2);
        b = b + (F(s+1)*(m*F(s+1)+b1)*sinc(2*k*F(s+1)) ...
                   - F(s)*(m*F(s)+b1)*sinc(2*k*F(s))) ...
                     * abs(W((s+1)/2)^2);
        if need_matrix
            G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))+sinc(2*I2*F(s+1))) ...
                   - .5*F(s)*(sinc(2*I1*F(s))+sinc(2*I2*F(s))) ) ...
                         * abs(W((s+1)/2)^2);
        end
    end;
    if Nodd
        b=[b0; b];
    end;

    if need_matrix
        a=G\b;
    else
        a=(W(1)^2)*4*b;
        if Nodd
            a(1) = a(1)/2;
        end
    end
    if Nodd
        h=[a(L+1:-1:2)/2; a(1); a(2:L+1)/2].';
    else
        h=.5*[flipud(a); a].';
    end;
elseif (ftype == 1),  % Type III and Type IV linear phase FIR
                  %  basis vectors are sin(2*pi*m*f) (see m below)
    if (differ),      % weight non-zero bands with 1/f^2
        do_weight = ( abs(M(1:2:length(M))) +  abs(M(2:2:length(M))) ) > 0;
    else
        do_weight = zeros(size(F));
    end

    if Nodd
        m=(1:L);      % type III
    else
        m=(0:L)+.5;   % type IV
    end;
    k=m';
    b=zeros(size(k));

    need_matrix = (~fullband) | (any(do_weight)) | (~constant_weights); 
    if need_matrix
        I1=k(:,ones(size(m)))+m(ones(size(k)),:);    % entries are m + k
        I2=k(:,ones(size(m)))-m(ones(size(k)),:);    % entries are m - k
        G=zeros(size(I1));
    end

    i = sqrt(-1);
    for s=1:2:length(F),
        if (do_weight((s+1)/2)),      % weight bands with 1/f^2
            if F(s) == 0, F(s) = 1e-5; end     % avoid singularities
            m=(M(s+1)-M(s))/(F(s+1)-F(s));
            b1=M(s)-m*F(s);
            snint1 = sineint(2*pi*k*F(s+1)) - sineint(2*pi*k*F(s));
            %snint1 = (-1/2/i)*(expint(i*2*pi*k*F(s+1)) ...
            %    -expint(-i*2*pi*k*F(s+1)) -expint(i*2*pi*k*F(s)) ...
            %    +expint(-i*2*pi*k*F(s)) );
            % csint1 = cosint(2*pi*k*F(s+1)) - cosint(2*pi*k*F(s)) ;
            csint1 = (-1/2)*(expint(i*2*pi*k*F(s+1))+expint(-i*2*pi*k*F(s+1))...
                            -expint(i*2*pi*k*F(s))  -expint(-i*2*pi*k*F(s)) );
            b=b + ( m*snint1 ...
                + b1*2*pi*k.*( -sinc(2*k*F(s+1)) + sinc(2*k*F(s)) + csint1 ))...
                * abs(W((s+1)/2)^2);
            snint1 = sineint(2*pi*F(s+1)*(-I2));
            snint2 = sineint(2*pi*F(s+1)*I1);
            snint3 = sineint(2*pi*F(s)*(-I2));
            snint4 = sineint(2*pi*F(s)*I1);
            G = G - ( ( -1/2*( cos(2*pi*F(s+1)*(-I2))/F(s+1)  ...
                      - 2*snint1*pi.*I2 ...
                      - cos(2*pi*F(s+1)*I1)/F(s+1) ...
                      - 2*snint2*pi.*I1 )) ...
                  - ( -1/2*( cos(2*pi*F(s)*(-I2))/F(s)  ...
                      - 2*snint3*pi.*I2 ...
                      - cos(2*pi*F(s)*I1)/F(s) ...
                      - 2*snint4*pi.*I1) ) ) ...
                         * abs(W((s+1)/2)^2);
        else      % use usual weights
            m=(M(s+1)-M(s))/(F(s+1)-F(s));
            b1=M(s)-m*F(s);
            b=b+(m/(4*pi*pi)*(sin(2*pi*k*F(s+1))-sin(2*pi*k*F(s)))./(k.*k))...
                         * abs(W((s+1)/2)^2) ;
            b = b + (((m*F(s)+b1)*cos(2*pi*k*F(s)) - ...
                 (m*F(s+1)+b1)*cos(2*pi*k*F(s+1)))./(2*pi*k)) ...
                         * abs(W((s+1)/2)^2) ;
            if need_matrix
                G = G + (.5*F(s+1)*(sinc(2*I1*F(s+1))-sinc(2*I2*F(s+1))) ...
                   - .5*F(s)*(sinc(2*I1*F(s))-sinc(2*I2*F(s)))) * ...
                      abs(W((s+1)/2)^2);
            end
        end;
    end

    if need_matrix
        a=G\b;
    else
        a=-4*b*(W(1)^2);
    end
    if Nodd
        h=.5*[flipud(a); 0; -a].';
    else
        h=.5*[flipud(a); -a].';
    end
    if differ, h=-h; end
end    

if nargout > 1
    a = 1;
end

function y = sineint(x)
% SINEINT (a.k.a. SININT)   Numerical Sine Integral
%   Used by FIRLS in the Signal Processing Toolbox.
%   Untested for complex or imaginary inputs.
%
%   See also SININT in the Symbolic Toolbox.

%   Was Revision: 1.5, Date: 1996/03/15 20:55:51

    i1 = find(real(x)<0);   % this equation is not valid if x is in the
                            % left-hand plane of the complex plane.
            % use relation Si(-z) = -Si(z) in this case (Eq 5.2.19, Abramowitz
            %  & Stegun).
    x(i1) = -x(i1);
    y = zeros(size(x));
    ind = find(x);
    % equation 5.2.21 Abramowitz & Stegun
    %  y(ind) = (1/(2*i))*(expint(i*x(ind)) - expint(-i*x(ind))) + pi/2;
    y(ind) = imag(expint(i*x(ind))) + pi/2;
    y(i1) = -y(i1);


function y=sinc(x)
%SINC Sin(pi*x)/(pi*x) function.
%   SINC(X) returns a matrix whose elements are the sinc of the elements 
%   of X, i.e.
%        y = sin(pi*x)/(pi*x)    if x ~= 0
%          = 1                   if x == 0
%   where x is an element of the input matrix and y is the resultant
%   output element.  

%   Author(s): T. Krauss, 1-14-93
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%       $Revision: 1.9 $  $Date: 1997/11/26 20:13:12 $

y=ones(size(x));
i=find(x);
y(i)=sin(pi*x(i))./(pi*x(i));
